setwd("/Users/fernandolucasruiz/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/TFM")## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
## Warning: package 'Hmisc' was built under R version 4.3.2
## Warning: package 'AnnotationDbi' was built under R version 4.3.1
## Warning: package 'IRanges' was built under R version 4.3.1
## Warning: package 'S4Vectors' was built under R version 4.3.1
library(org.Hs.eg.db)
library(VennDiagram)
library(ggplot2)
library(gridExtra)
library(purrr)
library(patchwork)
library(nortest)
library(ggrepel)
library(Rtsne)
library(umap)Descargar los datos
# covariables
ROSMAP_RINPMIAGESEX_covs <- readRDS("~/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/TFM/datos/ROSMAP_RINPMIAGESEX_covs.rds")
covs <- ROSMAP_RINPMIAGESEX_covs
#datos corregidos
ROSMAP_RINPMIAGESEX_resids <- readRDS("~/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/TFM/datos/ROSMAP_RINPMIAGESEX_resids.rds")
mat_exp <- ROSMAP_RINPMIAGESEX_residscovs$projid <- as.character(covs$projid)
covs$study <- as.factor(covs$study)
covs$ceradsc <- as.factor(covs$ceradsc)
covs$cogdx <- as.factor(covs$cogdx)
covs$neuroStatus <- as.factor(covs$neuroStatus)
rownames(covs) <- covs$mrna_id
describe(covs)## covs
##
## 13 Variables 342 Observations
## --------------------------------------------------------------------------------
## mrna_id
## n missing distinct
## 342 0 342
##
## lowest : 02_120405_2 04_120405_1 05_120405_4 07_120410_4 08_120410_5
## highest: 950_131107_8 951_131107_8 952_131107_8 954_131107_7 958_131107_7
## --------------------------------------------------------------------------------
## projid
## n missing distinct
## 342 0 342
##
## lowest : 10100574 10101327 10101589 10208143 10222853
## highest: 94722642 94828877 95491648 9841821 98953007
## --------------------------------------------------------------------------------
## study
## n missing distinct
## 342 0 2
##
## Value MAP ROS
## Frequency 167 175
## Proportion 0.488 0.512
## --------------------------------------------------------------------------------
## age_death
## n missing distinct Info Mean Gmd .05 .10
## 342 0 20 0.902 86.52 4.746 76.05 79.00
## .25 .50 .75 .90 .95
## 85.00 89.00 90.00 90.00 90.00
##
## Value 67 72 73 74 75 76 77 78 79 80 81
## Frequency 1 6 3 1 2 5 10 5 6 5 10
## Proportion 0.003 0.018 0.009 0.003 0.006 0.015 0.029 0.015 0.018 0.015 0.029
##
## Value 82 83 84 85 86 87 88 89 90
## Frequency 8 11 10 20 16 22 20 24 157
## Proportion 0.023 0.032 0.029 0.058 0.047 0.064 0.058 0.070 0.459
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## msex
## n missing distinct
## 342 0 2
##
## Value 0 1
## Frequency 223 119
## Proportion 0.652 0.348
## --------------------------------------------------------------------------------
## batch
## n missing distinct
## 342 0 9
##
## Value 0 1 2 3 4 5 6 7 8
## Frequency 11 55 48 47 44 54 38 30 15
## Proportion 0.032 0.161 0.140 0.137 0.129 0.158 0.111 0.088 0.044
## --------------------------------------------------------------------------------
## pmi
## n missing distinct Info Mean Gmd .05 .10
## 342 0 196 1 7.378 4.496 3.000 3.557
## .25 .50 .75 .90 .95
## 4.454 6.000 8.562 13.663 18.167
##
## lowest : 1 1.33333 1.75 1.83333 1.95
## highest: 21.4167 22.25 23.3333 23.9167 26
## --------------------------------------------------------------------------------
## rin
## n missing distinct Info Mean Gmd .05 .10
## 342 0 42 0.999 7.049 1.149 5.300 5.600
## .25 .50 .75 .90 .95
## 6.300 7.200 7.875 8.300 8.500
##
## lowest : 5 5.1 5.2 5.3 5.4, highest: 8.7 8.8 9 9.1 9.2
## --------------------------------------------------------------------------------
## braaksc
## n missing distinct
## 342 0 7
##
## Value 0 1 2 3 4 5 6
## Frequency 6 34 29 76 94 96 7
## Proportion 0.018 0.099 0.085 0.222 0.275 0.281 0.020
## --------------------------------------------------------------------------------
## ceradsc
## n missing distinct
## 342 0 4
##
## Value 1 2 3 4
## Frequency 108 64 38 132
## Proportion 0.316 0.187 0.111 0.386
## --------------------------------------------------------------------------------
## cogdx
## n missing distinct
## 342 0 5
##
## Value 1 2 3 4 5
## Frequency 102 62 4 159 15
## Proportion 0.298 0.181 0.012 0.465 0.044
## --------------------------------------------------------------------------------
## apoe_genotype
## n missing distinct
## 342 0 6
##
## Value 22 23 24 33 34 44
## Frequency 3 37 7 208 83 4
## Proportion 0.009 0.108 0.020 0.608 0.243 0.012
## --------------------------------------------------------------------------------
## neuroStatus
## n missing distinct
## 342 0 2
##
## Value 0 1
## Frequency 168 174
## Proportion 0.491 0.509
## --------------------------------------------------------------------------------
# geneset de Alzheimer extraidos de KEGG
# library(limma)
# library(AnnotationDbi)
# library(org.Hs.eg.db)
tab <- getGeneKEGGLinks(species="hsa")
tab$Symbol <- mapIds(org.Hs.eg.db, tab$GeneID,
column="SYMBOL", keytype="ENTREZID")## 'select()' returned 1:1 mapping between keys and columns
Vamos a seleccionar solamente los genes de Alzheimer de KEGG en nuestro dataset de expresión génica.
Visualizamos si se han seleccionado todos los genes y cuantos del geneset de KEGG
#library(VennDiagram)
#library(grid)
venn.plot <- venn.diagram(
x = list(GenesMatriz = colnames(mat_exp), GenesetAlzheimer = geneset_alz),
category.names = c("Matrix Genes", "Geneset Alzheimer"),
filename = NULL,
output = FALSE, # Asegura que no se exporte a un archivo
fill = c("#00CED1", "#E9967A"),
cex = 1, # Aumenta el tamaño del texto
fontface = "bold",
cat.cex = 1, # Aumenta el tamaño del texto de las categorÃas
cat.fontface = "bold",
cat.default.pos = "text",
cat.pos = 25, #posicion de las categoricas
cat.dist = 0.1, #distancia de las categoricas
rotation.degree = 0,
margin = 0.1, # hacerla más pequeña
lwd=0.5,
lty = "dashed", # Estilo de lÃnea discontinua
edge.col = "grey", # Color de los bordes
main = "Alzheimer genes in exp matrix",
main.fontface= "bold",
main.cex = 2,
main.pos = c(0.5, 1)
)
grid.newpage() # Asegura que el lienzo esté limpio
grid.draw(venn.plot)Hacemos PCA de nuestra nueva matriz con los genes seleccionados de Alzheimer. Los datos vienen libres de sesgos. Tampoco serÃa necesario el escalado de los datos ya que todos tienen la misma escala y vienen normalizados de antes.
El 63% de la varinza se explica en la primera componente principal. Eso en una matriz de expresión es raro, ya que suele ser bastante homogéneo.
pca.covs <- as.data.frame(pca_alz_genes$x[,1:10])
pca.covs <- merge(pca.covs, as.data.frame(covs), by = "row.names")
rownames(pca.covs) <- pca.covs$Row.names
pca.covs$Row.names <- NULL
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 11:ncol(pca.covs)) {
# Crea una gráfica para cada variable
p <- ggplot(pca.covs, aes_string(x = names(pca.covs)[1], y = names(pca.covs)[2], color = names(pca.covs)[i])) +
geom_point() +
labs(title = names(pca.covs)[i],
x = names(pca.covs)[1], y = names(pca.covs)[2], color = names(pca.covs)[i]) +
modelr::geom_ref_line(h = 0) +
modelr::geom_ref_line(v = 0) +
geom_point() +
xlab("PC 1 ") +
ylab("PC 2")
plots[[colnames(pca.covs)[i]]] <- p
}
grid.arrange(plots$study, plots$age_death, plots$msex, plots$batch, plots$pmi, plots$rin, ncol=2)grid.arrange(plots$braaksc, plots$ceradsc, plots$cogdx, plots$apoe_genotype, plots$neuroStatus, ncol=2)
### Seleccionar muestras con el doble de la desviación tÃpica
Seleccionamos muestras que esten 2 veces la desviacion tipica para la PC1 y PC2
outlierspc1 <- as.data.frame(pca_alz_genes$x[abs(pca_alz_genes$x[,1]) > 2*(pca_alz_genes$sdev[1]),1])
outlierspc2 <- as.data.frame(pca_alz_genes$x[abs(pca_alz_genes$x[,2]) > 2*(pca_alz_genes$sdev[2]),2])Las ploteamos
Creamos los grupos de muestras seleccionados de los outliers de PC1 y/o PC2
mPC1 <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc1),]
mPC2 <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc2),]
mPC12 <- mat_exp[(rownames(mat_exp) %in% rownames(outlierspc2)) & (rownames(mat_exp) %in% rownames(outlierspc1)),]
not_in_PC12 <- mat_exp[!((rownames(mat_exp) %in% rownames(outlierspc2)) | (rownames(mat_exp) %in% rownames(outlierspc1))),]# Asegurarse de que todos son vectores del mismo largo para el dataframe final
max_length <- max(length(rownames(mPC1)), length(rownames(mPC2)), length(rownames(mPC12)))
# Normalizar la longitud de los vectores (en caso de que alguno sea más corto)
rownames_mPC1 <- c(rownames(mPC1), rep("", max_length - length(rownames(mPC1))))
rownames_mPC2 <- c(rownames(mPC2), rep("", max_length - length(rownames(mPC2))))
rownames_mPC12 <- c(rownames(mPC12), rep("", max_length - length(rownames(mPC12))))
final_table <- data.frame(
mPC1 = rownames_mPC1,
mPC2 = rownames_mPC2,
mPC1_and_PC2 = rownames_mPC12
)
print(final_table)## mPC1 mPC2 mPC1_and_PC2
## 1 500_120515_2 681_120604_1 500_120515_2
## 2 629_120524_2 380_120503_1 629_120524_2
## 3 507_120515_2 500_120515_2 507_120515_2
## 4 605_120523_2 405_120503_2 271_120430_4
## 5 271_120430_4 629_120524_2 367_120502_4
## 6 367_120502_4 507_120515_2
## 7 657_120529_3
## 8 316_120501_3
## 9 584_120522_4
## 10 271_120430_4
## 11 367_120502_4
## 12 337_120501_5
## 13 232_120425_5
## 14 234_120425_6
## 15 790_130530_7
venn.plot <- venn.diagram(
x = list(mat_exp = rownames(pca_alz_genes$x), mPC1 = rownames(mPC1), mPC2 = rownames(mPC2)),
category.names = c("Expression matrix", "mPC1", "mPC2" ),
filename = NULL,
output = FALSE, # Asegura que no se exporte a un archivo
fill = c("#F08080", "#7FFFD4", "#87CEEB"),
cex = 1, # Aumenta el tamaño del texto
fontface = "bold",
cat.cex = 1, # Aumenta el tamaño del texto de las categorÃas
cat.fontface = "bold",
cat.default.pos = "text",
cat.pos = 0, #posicion de las categoricas
cat.dist = -0.05, #distancia de las categoricas
rotation.degree = 0,
margin = 0.1, # hacerla más pequeña
lwd=0.5,
lty = "dashed", # Estilo de lÃnea discontinua
edge.col = "grey", # Color de los bordes
main = "Sample set not scaled PCA",
main.fontface= "bold",
main.cex = 2,
main.pos = c(0.5, 1)
)
grid.newpage() # Asegura que el lienzo esté limpio
grid.draw(venn.plot)Metemos una nueva variable en las covariables dependiendo de si son outliers de PC1 o del PC2 o de ninguno.
rownames(covs) <- covs$mrna_id
for (i in rownames(mat_exp)){
if (i %in% rownames(mPC1) & !(i %in% rownames(mPC12))){
covs[i, "sampleset_PCA"] <- "mPC1"
}
if (i %in% rownames(mPC2) & !(i %in% rownames(mPC12))){
covs[i, "sampleset_PCA"] <- "mPC2"
}
if (i %in% rownames(mPC12)){
covs[i, "sampleset_PCA"] <- "mPC1 and mPC2"
}
if (!(i %in% rownames(mPC1)) & !(i %in% rownames(mPC2))){
covs[i, "sampleset_PCA"] <- "Not in both"
}
}
covs$sampleset_PCA <- as.factor(covs$sampleset_PCA)
table(covs$sampleset_PCA)##
## mPC1 mPC1 and mPC2 mPC2 Not in both
## 1 5 10 326
Ahora vamos a hacer un análisis estadÃstico para ver si hay diferencias significativas entre las muestras outliers de cada componente principal con las que no lo están (mPC1 vs Not in both, mPC2 vs Not in both, mPC12 vs Not in both)
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 1:ncol(covs)) {
col_name <- colnames(covs)[i]
if (class(covs[[i]]) == "character"){
next
}
if (class(covs[[i]]) == "factor"){
p <- covs %>%
group_by(across(all_of(col_name)), sampleset_PCA) %>%
summarise(count = n(), .groups = 'drop') %>%
mutate(total = sum(count), percentage = (count / total) * 100) %>%
ggplot(aes_string(x = names(covs)[i], y = "percentage", fill = "sampleset_PCA")) +
geom_bar(stat = "identity", position = "fill") +
labs(x = NULL, y = "%", fill = col_name) +
theme_minimal()
plots[[colnames(covs)[i]]] <- p
}
if (class(covs[[i]]) == "numeric"){
p <- ggplot(covs, aes_string(y=names(covs)[i], x = "sampleset_PCA", fill = "sampleset_PCA")) +
geom_boxplot(alpha=.7)+
geom_jitter(color ="#8B1A1A", size = 0.7, alpha = 0.7) +
theme_minimal() +
labs(title = names(covs)[i], x = names(covs)[i])
plots[[colnames(covs)[i]]] <- p
}
}
grid.arrange(plots$age_death, plots$pmi, plots$rin, ncol=2)statistical.tests <- function(df, outlier){
for (i in colnames(df)) {
#outlier <- "mPC1"
if (i == "mrna_id" | i == "projid" | i == "sampleset_PCA" | i == "sampleset_tSNE") {
next
}
# filtro primero los datos por los que quiero comparar outlier y no incluidos
filtered_df <- df[df$sampleset_PCA %in% c(outlier, "Not in both"), ]
# si la variable es numérica
if (class(filtered_df[[i]]) == "numeric") {
# hacemos test de normalidad de la variable para ver si hacemos test paramétrico o no paramétrico
normality <- lillie.test(filtered_df[[i]])
# si los datos no son paramétricos
if (normality$p.value < 0.05) {
# Verificar que sampleset_PCA tenga exactamente 2 niveles después del filtrado
if (length(unique(filtered_df$sampleset_PCA)) == 2) {
wilcox_test <- wilcox.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA)
if (wilcox_test$p.value < 0.05) {
print(paste("U Mann-Witney test significativo de", i, "con un pvalor de:", round(wilcox_test$p.value, 6), "entre", outlier, " y 'Not in both'"))
} else {
print(paste("No hay diferencia significativa de la variable", i, "entre", outlier, "y 'Not in both'"))
}
}
# si la variable numérica es paramétrica
} else if (normality$p.value >= 0.05) {
# hacemos test de la varianza
var_test <- var.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA)
# le pasamos el t.test dependiendo de si las varanzas son iguales o no
if (var_test$p.value < 0.05) {
t_test <- t.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA, var.equal = F)
if (t_test$p.value < 0.05){
print(paste("T test significativo (varianzas desiguales) de",
i, "con un pvalor de:", round(t_test$p.value, 6)))
} else {
print(paste("No hay diferencia significativa de la variable", i, "entre", outlier, "y 'Not in both'"))
}
} else {
t_test <- t.test(filtered_df[[i]] ~ filtered_df$sampleset_PCA, var.equal = T)
if (t_test$p.value < 0.05){
print(paste("T test significativo (varianzas iguales) de",
i, "con un pvalor de:", round(t_test$p.value, 6)))
} else {
print(paste("No hay diferencia significativa de la variable", i, "entre", outlier, "y 'Not in both'"))
}
}
}
}
# si la variable es categórica le vamos a pasar un test chisq
if (class(filtered_df[[i]]) == "factor") {
# hacemos la tabla de contingencia
contingency_table <- table(filtered_df[[i]], filtered_df$sampleset_PCA)
# test de Chi-cuadrado
test_chi <- chisq.test(contingency_table)
if (!is.na(test_chi$p.value) && test_chi$p.value < 0.05) {
print(paste("Test de Chi-cuadrado significativo para", i, "con un p-valor de:", round(test_chi$p.value, 6), "entre", outlier, " y 'Not in both'"))
} else {
print(paste("No hay diferencia significativa de la variable", i, "entre", outlier, "y 'Not in both'"))
}
}
}
}## [1] "No hay diferencia significativa de la variable study entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable study entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable study entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC12 y 'Not in both'"
# Scree plot con los datos escalados
var_exp <- pca_alz_genes_scaled$sdev^2
prop_var_exp <- var_exp / sum(var_exp)
cum_var_exp <- cumsum(prop_var_exp)
df_var_exp <- data.frame(Comp = 1:length(prop_var_exp), VarExp = prop_var_exp)
df_cum_var_exp <- data.frame(Comp = 1:length(cum_var_exp), CumVarExp = cum_var_exp)
ggplot(df_var_exp[1:20,], aes(x = Comp, y = VarExp)) +
geom_bar(stat = "identity", fill = "skyblue") +
geom_line(aes(group = 1), color = "blue") +
geom_point(color = "blue") +
theme_minimal() +
labs(x = "Principal components", y = "Variance", title = "Scree Plot scaled") +
ylim(c(0,1)) +
geom_line(data = df_cum_var_exp[1:20,], aes(x = Comp, y = CumVarExp), color="#8B1A1A") +
geom_point(data = df_cum_var_exp[1:20,], aes(x = Comp, y = CumVarExp), color = "red") +
geom_bar(data = df_cum_var_exp[1:20,], aes(x = Comp, y = CumVarExp), stat = "identity", fill = "red", alpha= 0.25) +
annotate("text", x = 4, y = 0.85, label = "Cumulative Scree Plot", color = "#8B1A1A", size = 4) +
geom_text(data = df_cum_var_exp[1:20,], aes(x=Comp, y = CumVarExp +0.04, label = round(CumVarExp, 2)))
Vemos que con las 20 primeras no encontramos ni siquiera el 80% de la
varianza.
pca.scaled.covs <- as.data.frame(pca_alz_genes_scaled$x[,1:10])
pca.scaled.covs <- merge(pca.scaled.covs, as.data.frame(covs), by = "row.names")
rownames(pca.scaled.covs) <- pca.scaled.covs$Row.names
pca.scaled.covs$Row.names <- NULL
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 11:ncol(pca.scaled.covs)) {
# Crea una gráfica para cada variable
p <- ggplot(pca.scaled.covs, aes_string(x = names(pca.scaled.covs)[1], y = names(pca.scaled.covs)[2], color = names(pca.scaled.covs)[i])) +
geom_point() +
labs(title = names(pca.scaled.covs)[i],
x = names(pca.scaled.covs)[1], y = names(pca.scaled.covs)[2], color = names(pca.scaled.covs)[i]) +
modelr::geom_ref_line(h = 0) +
modelr::geom_ref_line(v = 0) +
geom_point() +
xlab("PC 1 ") +
ylab("PC 2")
plots[[colnames(pca.scaled.covs)[i]]] <- p
}
grid.arrange(plots$study, plots$age_death, plots$msex, plots$batch, plots$pmi, plots$rin, ncol=2)grid.arrange(plots$braaksc, plots$ceradsc, plots$cogdx, plots$apoe_genotype, plots$neuroStatus, ncol=2)outlierspc1_scaled <- as.data.frame(pca_alz_genes_scaled$x[abs(pca_alz_genes_scaled$x[,1]) > 2*(pca_alz_genes_scaled$sdev[1]),1])
outlierspc2_scaled <- as.data.frame(pca_alz_genes_scaled$x[abs(pca_alz_genes_scaled$x[,2]) > 2*(pca_alz_genes_scaled$sdev[2]),2])df <- as.data.frame(pca_alz_genes_scaled$x)
plot1 <- ggplot(df, aes(x = PC1)) +
geom_density(fill = "#00CED1", alpha = 0.5) +
geom_vline(xintercept = mean(df$PC1) + 2*pca_alz_genes_scaled$sdev[1], linetype = "dashed", color = "blue") +
geom_vline(xintercept = mean(df$PC1) - 2*pca_alz_genes_scaled$sdev[1], linetype = "dashed", color = "blue") +
labs(title = "PC1 density with 2*SD scaled PCA") +
geom_text(data = outlierspc1_scaled,
aes(x = outlierspc1_scaled[,1], y= 0, label = rownames(outlierspc1_scaled)),
vjust = 1.5, hjust = 0, size = 3, color = "#E9965A", angle = 90, fontface= "bold") +
geom_rug(data = as.data.frame(pca_alz_genes_scaled$x), aes(x= PC1, y = 0), color= ifelse(abs(pca_alz_genes_scaled$x[,1]) > 2*(pca_alz_genes_scaled$sdev[1]), "#8B1A1A", "grey")) +
theme_minimal()
plot2 <- ggplot(df, aes(x = PC2)) +
geom_density(fill = "#00CED1", alpha = 0.5) +
geom_vline(xintercept = mean(df$PC2) + 2*pca_alz_genes_scaled$sdev[2], linetype = "dashed", color = "blue") +
geom_vline(xintercept = mean(df$PC2) - 2*pca_alz_genes_scaled$sdev[2], linetype = "dashed", color = "blue") +
labs(title = "PC2 density with 2*SD scaled PCA") +
geom_text(data = outlierspc2_scaled,
aes(x = outlierspc2_scaled[,1], y= 0, label = rownames(outlierspc2_scaled)),
vjust = 1.5, hjust = 0, size = 3, color = "#E9965A", angle = 90, fontface= "bold") +
geom_rug(data = as.data.frame(pca_alz_genes_scaled$x), aes(x= PC2, y = 0), color= ifelse(abs(pca_alz_genes_scaled$x[,2]) > 2*(pca_alz_genes_scaled$sdev[2]), "#8B1A1A", "grey")) +
theme_minimal()
grid.arrange(plot1, plot2, ncol = 1)mPC1_scaled <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc1_scaled),]
mPC2_scaled <- mat_exp[rownames(mat_exp) %in% rownames(outlierspc2_scaled),]
mPC12_scaled <- mat_exp[(rownames(mat_exp) %in% rownames(outlierspc2_scaled)) & (rownames(mat_exp) %in% rownames(outlierspc1_scaled)),]
not_in_PC12_scaled <- mat_exp[!((rownames(mat_exp) %in% rownames(outlierspc2_scaled)) | (rownames(mat_exp) %in% rownames(outlierspc1_scaled))),]# Asegurarse de que todos son vectores del mismo largo para el dataframe final
max_length <- max(length(rownames(mPC1_scaled)), length(rownames(mPC2_scaled)), length(rownames(mPC12_scaled)))
# Normalizar la longitud de los vectores (en caso de que alguno sea más corto)
rownames_mPC1_scaled <- c(rownames(mPC1_scaled), rep("", max_length - length(rownames(mPC1_scaled))))
rownames_mPC2_scaled <- c(rownames(mPC2_scaled), rep("", max_length - length(rownames(mPC2_scaled))))
rownames_mPC12_scaled <- c(rownames(mPC12_scaled), rep("", max_length - length(rownames(mPC12_scaled))))
final_table <- data.frame(
mPC1 = rownames_mPC1_scaled,
mPC2 = rownames_mPC2_scaled,
mPC1_and_PC2 = rownames_mPC12_scaled
)
print(final_table)## mPC1 mPC2 mPC1_and_PC2
## 1 274_120430_1 384_120503_1
## 2 724_120605_2 380_120503_1
## 3 534_120516_2 300_120430_2
## 4 657_120529_3 07_120410_4
## 5 193_120424_3 186_120424_4
## 6 584_120522_4 628_120524_4
## 7 271_120430_4 336_120501_5
## 8 337_120501_5 691_120605_5
## 9 79_120417_5 08_120410_5
## 10 215_120425_5 119_120418_6
## 11 232_120425_5 893_130923_7
## 12 505_120515_6 905_131010_7
## 13 537_120516_6
## 14 769_130523_7
## 15 929_131031_8
venn.plot <- venn.diagram(
x = list(mat_exp = rownames(mat_exp), mPC1 = rownames(mPC1_scaled), mPC2 = rownames(mPC2_scaled)),
category.names = c("Expression matrix", "mPC1", "mPC2" ),
filename = NULL,
output = FALSE, # Asegura que no se exporte a un archivo
fill = c("#F08080", "#7FFFD4", "#87CEEB"),
cex = 1, # Aumenta el tamaño del texto
fontface = "bold",
cat.cex = 1, # Aumenta el tamaño del texto de las categorÃas
cat.fontface = "bold",
cat.default.pos = "text",
cat.pos = 0, #posicion de las categoricas
cat.dist = -0.05, #distancia de las categoricas
rotation.degree = 0,
margin = 0.1, # hacerla más pequeña
lwd=0.5,
lty = "dashed", # Estilo de lÃnea discontinua
edge.col = "grey", # Color de los bordes
main = "Sample set scaled PCA",
main.fontface= "bold",
main.cex = 2,
main.pos = c(0.5, 1)
)
grid.newpage() # Asegura que el lienzo esté limpio
grid.draw(venn.plot)Voy ahora a ver de esas muestras seleccionadas, qué covariables están asociadas con esas muestras outliers
rownames(covs) <- covs$mrna_id
covs.scaled <- covs
for (i in rownames(covs.scaled)){
if (i %in% rownames(mPC1_scaled) & !(i %in% rownames(mPC12_scaled))){
covs.scaled[i, "sampleset_PCA"] <- "mPC1"
}
if (i %in% rownames(mPC2_scaled) & !(i %in% rownames(mPC12_scaled))){
covs.scaled[i, "sampleset_PCA"] <- "mPC2"
}
if (i %in% rownames(mPC12_scaled)){
covs.scaled[i, "sampleset_PCA"] <- "mPC1 and mPC2"
}
if (!(i %in% rownames(mPC1_scaled)) & !(i %in% rownames(mPC2_scaled))){
covs.scaled[i, "sampleset_PCA"] <- "Not in both"
}
}
covs.scaled$sampleset_PCA <- as.factor(covs.scaled$sampleset_PCA)
table(covs.scaled$sampleset_PCA)##
## mPC1 mPC1 and mPC2 mPC2 Not in both
## 15 0 12 315
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 1:ncol(covs.scaled)) {
col_name <- colnames(covs.scaled)[i]
if (class(covs.scaled[[i]]) == "character"){
next
}
if (class(covs.scaled[[i]]) == "factor"){
p <- covs.scaled %>%
group_by(across(all_of(col_name)), sampleset_PCA) %>%
summarise(count = n(), .groups = 'drop') %>%
mutate(total = sum(count), percentage = (count / total) * 100) %>%
ggplot(aes_string(x = names(covs.scaled)[i], y = "percentage", fill = "sampleset_PCA")) +
geom_bar(stat = "identity", position = "fill") +
labs(x = NULL, y = "%", fill = col_name) +
theme_minimal()
plots[[colnames(covs.scaled)[i]]] <- p
}
if (class(covs.scaled[[i]]) == "numeric"){
p <- ggplot(covs.scaled, aes_string(y=names(covs.scaled)[i], x = "sampleset_PCA", fill = "sampleset_PCA")) +
geom_boxplot(alpha=.7)+
geom_jitter(color ="#8B1A1A", size = 0.7, alpha = 0.7) +
theme_minimal() +
labs(title = names(covs.scaled)[i], x = names(covs.scaled)[i])
plots[[colnames(covs.scaled)[i]]] <- p
}
}
grid.arrange(plots$age_death, plots$pmi, plots$rin, ncol=2)## [1] "No hay diferencia significativa de la variable study entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable study entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable age_death entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable pmi entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable rin entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable study entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mPC12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mPC12 y 'Not in both'"
Parámetros clave en t-SNE: Perplexity (Perplejidad): Este es uno de los parámetros más importantes en t-SNE. La perplejidad se relaciona con el número de vecinos cercanos que t-SNE considera cuando mapea los datos. Un valor demasiado bajo hace que el modelo se concentre demasiado en la estructura local, mientras que un valor demasiado alto puede llevar a una visión global que pierde detalles. La elección óptima depende del tamaño del dataset, pero valores tÃpicos están entre 5 y 50.
Número de iteraciones: t-SNE es un algoritmo iterativo, y el número de iteraciones puede afectar la estabilidad de los resultados. Un número insuficiente de iteraciones puede resultar en una organización incompleta, mientras que demasiadas iteraciones pueden no proporcionar beneficios adicionales y aumentar el tiempo de cálculo.
Tasa de aprendizaje: Este parámetro controla el tamaño del paso en cada actualización durante la optimización. Una tasa de aprendizaje muy baja puede hacer que el algoritmo tarde mucho en converger, mientras que una tasa demasiado alta puede llevar a una convergencia pobre.
set.seed(1234)
tsne <- Rtsne(mat_exp_alz_genes, dims = 2, theta = 0.0, )
tsne.data <- as.data.frame(tsne$Y)
row.names(tsne.data) <- row.names(mat_exp_alz_genes)
tsne.data.covs <- merge(tsne.data, covs, by = "row.names")
tsne.data.covs$Row.names <- NULL
row.names(tsne.data.covs) <- tsne.data.covs$mrna_id
ggplot(tsne.data.covs, aes(x=V1, y=V2, color=neuroStatus)) +
geom_point(size=2) +
guides(colour=guide_legend(override.aes=list(size=6))) +
xlab("tSNE1") + ylab("tSNE2") +
ggtitle("t-SNE") +
theme_minimal() +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank()) +
scale_colour_brewer(palette = "Set1")rownames(tsne.data) <- rownames(mat_exp_alz_genes)
# Calcular la media y la desviación estándar para cada componente de t-SNE
mean.tsne1 <- mean(tsne.data[,1])
sd.tsne1 <- sd(tsne.data[,1])
mean.tsne2 <- mean(tsne.data[,2])
sd.tsne2 <- sd(tsne.data[,2])
# Identificar muestras a más de 2 desviaciones estándar de la media
outliers.tsne1 <- tsne.data[abs(tsne.data[,1] - mean.tsne1) > 2 * sd.tsne1, ]
outliers.tsne2 <-tsne.data[abs(tsne.data[,2] - mean.tsne2) > 2 * sd.tsne2, ]## Warning in mean.default(df$V2): argument is not numeric or logical: returning
## NA
## Warning in mean.default(df$V2): argument is not numeric or logical: returning
## NA
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_vline()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_vline()`).
mtSNE1 <- mat_exp[rownames(mat_exp) %in% rownames(outliers.tsne1),]
mtSNE2 <- mat_exp[rownames(mat_exp) %in% rownames(outliers.tsne2),]
mtSNE12 <- mat_exp[(rownames(mat_exp) %in% rownames(outliers.tsne1)) & (rownames(mat_exp) %in% rownames(outliers.tsne2)),]
not_in_tSNE12 <- mat_exp[!((rownames(mat_exp) %in% rownames(outliers.tsne1)) | (rownames(mat_exp) %in% rownames(outliers.tsne2))),]# Asegurarse de que todos son vectores del mismo largo para el dataframe final
max_length <- max(length(rownames(mtSNE1)), length(rownames(mtSNE2)), length(rownames(mtSNE12)))
# Normalizar la longitud de los vectores (en caso de que alguno sea más corto)
rownames_tsne1 <- c(rownames(mtSNE1), rep("", max_length - length(rownames(mtSNE1))))
rownames_tsne2 <- c(rownames(mtSNE2), rep("", max_length - length(rownames(mtSNE2))))
rownames_tsne12 <- c(rownames(mtSNE12), rep("", max_length - length(rownames(mtSNE12))))
final_table <- data.frame(
mtSNE1 = rownames_tsne1,
mtSNE2 = rownames_tsne2,
mtSNE1_and_mtSNE2 = rownames_tsne12
)
print(final_table)## mtSNE1 mtSNE2 mtSNE1_and_mtSNE2
## 1 681_120604_1 666_120530_1
## 2 790_130530_7 426_120507_1
## 3 430_120507_1
## 4 380_120503_1
## 5 500_120515_2
## 6 629_120524_2
## 7 507_120515_2
## 8 544_120516_2
## 9 203_120424_2
## 10 605_120523_2
## 11 568_120521_3
## 12 292_120430_3
## 13 186_120424_4
## 14 367_120502_4
## 15 327_120501_5
## 16 336_120501_5
## 17 257_120426_5
venn.plot <- venn.diagram(
x = list(mat_exp = rownames(tsne.data), mtSNE1 = rownames(mtSNE1), mtSNE2 = rownames(mtSNE2)),
category.names = c("Expression matrix", "mtSNE1", "mtSNE2" ),
filename = NULL,
output = FALSE, # Asegura que no se exporte a un archivo
fill = c("#F08080", "#7FFFD4", "#87CEEB"),
cex = 1, # Aumenta el tamaño del texto
fontface = "bold",
cat.cex = 1, # Aumenta el tamaño del texto de las categorÃas
cat.fontface = "bold",
cat.default.pos = "text",
cat.pos = 0, #posicion de las categoricas
cat.dist = -0.05, #distancia de las categoricas
rotation.degree = 0,
margin = 0.1, # hacerla más pequeña
lwd=0.5,
lty = "dashed", # Estilo de lÃnea discontinua
edge.col = "grey", # Color de los bordes
main = "Sample set tSNE",
main.fontface= "bold",
main.cex = 2,
main.pos = c(0.5, 1)
)
grid.newpage() # Asegura que el lienzo esté limpio
grid.draw(venn.plot)for (i in rownames(covs)){
if (i %in% rownames(mtSNE1) & !(i %in% rownames(mtSNE12))){
covs[i, "sampleset_tSNE"] <- "mtSNE1"
}
if (i %in% rownames(mtSNE2) & !(i %in% rownames(mtSNE12))){
covs[i, "sampleset_tSNE"] <- "mtSNE2"
}
if (i %in% rownames(mtSNE12)){
covs[i, "sampleset_tSNE"] <- "mtSNE1 and mtSNE2"
}
if (!(i %in% rownames(mtSNE1)) & !(i %in% rownames(mtSNE2))){
covs[i, "sampleset_tSNE"] <- "Not in both"
}
}
covs$sampleset_tSNE <- as.factor(covs$sampleset_tSNE)
table(covs$sampleset_tSNE)##
## mtSNE1 mtSNE2 Not in both
## 2 17 323
plots <- list()
# Itera sobre las columnas de covariables y añadir a la lista
for (i in 1:ncol(covs)) {
col_name <- colnames(covs)[i]
if (class(covs[[i]]) == "character"){
next
}
if (class(covs[[i]]) == "factor"){
p <- covs %>%
group_by(across(all_of(col_name)), sampleset_tSNE) %>%
summarise(count = n(), .groups = 'drop') %>%
mutate(total = sum(count), percentage = (count / total) * 100) %>%
ggplot(aes_string(x = names(covs)[i], y = "percentage", fill = "sampleset_tSNE")) +
geom_bar(stat = "identity", position = "fill") +
labs(x = NULL, y = "%", fill = col_name) +
theme_minimal()
plots[[colnames(covs)[i]]] <- p
}
if (class(covs[[i]]) == "numeric"){
p <- ggplot(covs, aes_string(y=names(covs)[i], x = "sampleset_tSNE", fill = "sampleset_tSNE")) +
geom_boxplot(alpha=.7)+
geom_jitter(color ="#8B1A1A", size = 0.7, alpha = 0.7) +
theme_minimal() +
labs(title = names(covs)[i], x = names(covs)[i])
plots[[colnames(covs)[i]]] <- p
}
}
grid.arrange(plots$age_death, plots$pmi, plots$rin, ncol=2)## [1] "No hay diferencia significativa de la variable study entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mtSNE1 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable study entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mtSNE2 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable study entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable msex entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable batch entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable braaksc entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable ceradsc entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable cogdx entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable apoe_genotype entre mtSNE12 y 'Not in both'"
## [1] "No hay diferencia significativa de la variable neuroStatus entre mtSNE12 y 'Not in both'"
local.config <- umap.defaults
# local.config$n_neighbors <- 4
# local.config$n_components <- 2
# local.config$n_epochs <- 100
# local.config$metric<- "euclidean"
set.seed(1234)
umap.ad <- umap(mat_exp_alz_genes,random_stage=1234, local.config)
umap.data <- as.data.frame(umap.ad$layout)
row.names(umap.data) <- row.names(mat_exp_alz_genes)
umap.data.covs <- merge(umap.data, covs, by = "row.names")
umap.data.covs$Row.names <- NULL
row.names(umap.data.covs) <- umap.data.covs$mrna_id
ggplot(umap.data.covs, aes(x=V1, y=V2, color = neuroStatus)) +
geom_point(size=2) +
guides(colour=guide_legend(override.aes=list(size=6))) +
xlab("UMAP1") + ylab("UMAP2") +
ggtitle("UMAP") +
theme_minimal() +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank()) +
scale_colour_brewer(palette = "Set2")